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Abstract 

We review recent results of molecular dynamics simulations of two 
models of liquid water, the extended simple point charge (SPC/E) and 
the Mahoney-Jorgensen transferable intermolecular potential with five 
points (TIP5P), which is closer to real water than previously-proposed 
classical pairwise additive potentials. Simulations of the TIP5P model 
for a wide range of deeply supercooled states, including both posi- 
tive and negative pressures, reveal (i) the existence of a non-monotonic 
"nose-shaped" temperature of maximum density (TMD) line and a non- 
reentrant spinodal, (ii) the presence of a low temperature phase tran- 
sition. The TMD that changes slope from negative to positive as P 
decreases and, notably, the point of crossover between the two behav- 
iors is located at ambient pressure (temperature ~ 4 C , and density 
~ 1 g/cm 3 ). Simulations on the dynamics of the SPC/E model reveal 

(iii) the dynamics at negative pressure shows a minimum in the diffu- 
sion constant D when the density is decreased at constant temperature, 
complementary to the known maximum of D at higher pressures, and 

(iv) the loci of minima of D relative to the spinodal shows that they are 
inside the thermodynamically metastable regions of the phase-diagram. 
These dynamical results reflect the initial enhancement and subsequent 
breakdown of the tetrahedral structure and of the hydrogen bond net- 
work as the density decreases. 
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1. Introduction 

Water is an important liquid in nature, and is also fundamental in 
chemical and technological applications. Although the individual water 
molecule has a simple chemical structure, water is considered a complex 
fluid because of its anomalous behavior [1, 2, 3, 4, 5, 6]. It expands 
on freezing and, at a pressure of 1 atm, the density has a maximum at 
4°C. Additionally, there is a minimum of the isothermal compressibil- 
ity at 46° C and a minimum of the isobaric heat capacity at 35° C [7]. 
These anomalies are linked with the microscopic structure of liquid wa- 
ter, which can be regarded as a transient gel — a highly associated liquid 
with strongly directional hydrogen bonds [8, 9]. Each water molecule 
acts as both a donor and an acceptor of bonds, generating a structure 
that is locally ordered, similar to that of ice, but maintaining the long- 
range disorder typical of liquids. Despite the extensive work that has 
been done on water, many aspects of its behavior remain unexplained. 

Several scenarios have been proposed to account for the the anoma- 
lous behavior of the thermodynamic response functions on cooling, each 
predicting a different behavior for the liquid spinodal, the line of the 
limit of stability separating the region where liquid water is metastable 
from the region where the liquid is unstable, (i) According the stability- 
limit conjecture [10, 11], the pressure of the spinodal line should decrease 
on cooling, become negative, and increase again after passing through a 
minimum. It reenters the positive pressure region of the phase diagram 
at a very low temperature, thereby giving rise to a line of singulari- 
ties in the positive pressure region, and consequently the increase in the 
thermodynamic response functions on cooling in the anomalous region 
is due to the proximity of this reentrant spinodal. (ii) The critical point 
hypothesis [12, 13, 14, 15, 16, 17], proposes a new critical point at the ter- 
minus of a first-order phase transition line separating two liquid phases 
of different density. The anomalous increases of the response functions, 
compressibility, specific heat, and volume expansivity, is interpreted in 
terms of this critical point, (iii) The singularity-free hypothesis [9, 18, 19] 
proposes that actually there is no divergence close to the anomalous re- 
gion; the response functions grow on lowering temperature but remain 
finite, attaining maximum values. 

Water properties and anomalies can be strongly influenced by the 
physical or chemical properties of the medium [1, 3, 4, 5, 20, 21]. The ef- 
fect not only of applied pressure, but also of negative pressure ( "stretch- 
ing" ) is remarkable. The study of the behavior of this fluid under nega- 
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tive pressures is relevant not only from the academic point of view, but 
also for realistic systems. For example, negative pressures are observed 
[22], and seem to play an important role in the mechanism of water 
transport in plants. Therefore, properties that modify the structure of 
water, especially if this modification is similar to the effect of stretch- 
ing (as is the case in some hydrogels [20]), also influence its dynamical 
behavior. 

Dynamic properties, such as the diffusion constant, have been studied 
in detail for water systems at atmospheric and at high positive pres- 
sures, both experimentally [23, 24] as well as by computer simulations 
[25, 26, 27, 28, 29, 30, 31]. The increase of pressure increases the pres- 
ence of defects and of interstitial water molecules in the network [26]. 
They disrupt the tetrahedral local structure, weakening the hydrogen 
bonds, and thus increasing the diffusion constant [30, 31]. However, a 
further increase in the pressure leads to steric effects which works in the 
direction of lowering the mobility. The interplay of these factors leads to 
a maximum in the diffusion constant [30, 31] at some high density /9 max . 
Above this density (or corresponding pressure) , the diffusion of water is 
in some sense like that of a normal liquid, controlled by hindrance, with 
the hydrogen bonds playing a secondary role. However, the behavior at 
very low p is less well understood. 

2. Location of the Spinodal at Positive and 
Negative Pressures 

Relatively few experimental works [32, 33] and simulations [12, 16, 
29, 30, 35, 34] have been performed on "stretched" water. In this nega- 
tive pressure region of the phase diagram the system is metastable, and 
becomes unstable beyond the spinodal line, so locating the spinodal we 
can ensure that the simulated state points lie in the metastable and not 
in the unstable region. Moreover, the shape of the spinodal can test the 
stability-limit conjecture against the critical point hypothesis and the 
singularity-free interpretation, so we first discuss the density and pres- 
sure of the spinodal, which we denote p sp (T) and P sp (T), respectively. 

Yamada and her coworkers [36] simulated a system of N = 343 
molecules interacting with the TIP5P potential [37]. TIP5P is a five- 
site, rigid, non-polarizable water model, not unlike the ST2 model [38]. 
The TIP5P potential accurately reproduces the density anomaly at 1 
atm and exhibits excellent structural properties when compared with 
experimental data [37, 39]. The TMD shows the correct pressure depen- 
dence, shifting to lower temperatures as pressure is increased. Under 
ambient conditions, the diffusion constant is close to the experimental 
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value, with reasonable temperature and pressure dependence away from 
ambient conditions [37]. Equilibration runs were performed at constant 
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Figure 1. Dependence on density of the pressure at all temperatures investigated 
(T = 215, 220, 230, 240, 250, 260, 270, 280, 290, 300, 320 K, from bottom to top). Each 
curve has been shifted by n x 150 MPa to avoid overlaps. An inflection appears as T 
is decreased, transforming into a "flat" coexistence region at T — 215 K, indicating 
the presence of a liquid-liquid transition. Inset: A detailed view of the T = 215 K 
isotherm. Adapted from [36]. 

T. After thermalization at T = 320 K the thermostat temperature 
was set to the temperature of interest. The system evolved for a time 
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longer than the structural relaxation time r a , defined as the time at 
which F s (Qo,T a ) = 1/e, where F s (Qo,t) is the self-intermediate scatter- 
ing function evaluated at Qq = 18 nm _1 , the location of the first peak 
of the static structure factor. In the time T a , each molecule diffuses on 
average a distance of the order of the nearest neighbor distance. We 
use the final configuration of the equilibration run to start a production 
run of length greater than several r Q and then analyze the calculated 
trajectory. 




Figure 2. Pressure along seven isochores; the minima correspond to the temperature 
of maximum density line (dashed line). Note the "nose" of the TMD line at T = 4 
C. Stars denote the liquid spinodal line, which is not reentrant, and terminates at the 
liquid-gas critical point. Adapted from [36]. 

Figure 1 shows results for pressure along isotherms. At lower tem- 
peratures an inflection develops, which becomes a "flat" isotherm at 
the lowest temperature, T = 215 K. The presence of a flat region in- 
dicates that a phase separation takes place; the critical temperature is 
T c > = 217 ± 3 K, the critical pressure is P c > = 340± 20 MPa, and the 
critical density pc = 1.13 ± 0.04 g/cm 3 . 

Figure 2(a) plots the pressure along isochores. The curves show min- 
ima as a function of temperature; the locus of the minima is the TMD 
line, since {dP/dT)y = ap/Kx- Note that the pressure exhibits a mini- 
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mum if the density passes through a maximum (a p = 0). It is clear that, 
as in the case of ST2 water, TIP5P water has a TMD that changes slope 
from negative to positive as P decreases. Notably, the point of crossover 
between the two behaviors is located at ambient pressure, ± ^ 4 C , and 
p ps 1 g/cm 3 . 

Also plotted the spinodal line, obtained by fitting the isotherms (for 
T > 300K) of Fig. 1 to the form P(T,p) = P S (T) + A [p - p s (T)] 2 , 
where P S (T) and p s {T) denote the pressure and density of the spinodal 
line. This functional form is the mean field prediction for P(p) close 
to a spinodal line. For T < 250K, P S (T) is calculated by estimating 
the location of the minimum of P(p). The results in Fig. 2 show that 
the liquid spinodal line is not reentrant and does not intersect the TMD 
line. 

3. Dynamic Properties 

We next discuss results on the dynamics of stretched water recently 
obtained by Netz and his collaborators [40]. While there are a large 
number of intermolecular potential functions used to simulate water, 
each of which gives slightly different results, the overall thermodynam- 
ics picture obtained from these models is generally very similar. Since 
dynamic properties are particularly sensitive to the potential choice, the 
extended simple point charge (SPC/E) potential is used since it repro- 
duces both the maximum in diffusivity under pressure as well as the 
power-law behavior of dynamics properties on cooling. For understand- 
ing the properties of water at negative pressure, simulations are partic- 
ularly important since experiments are very difficult to perform in this 
region. 

The effect of extreme conditions on the flow of the liquid is assessed 
by calculating the diffusion constant D, defined by the asymptotic value 
of the slope of the mean square displacement versus time. We show D 
along isotherms in Fig. 3. For T < 260 K, D has a minimum value at p ps 
0.9 g/cm 3 , which becomes more pronounced at lower T. This behavior 
can be understood considering the structural changes that occur with 
decreasing density. At low T, the decreased density enhances the local 
tetrahedral ordering, which leads to a decrease in D. Further decreases 
in density reduces the stability of the tetrahedral structure and causes 
an increase of D. 

The location of the minimum is near the ice Ih density ps 0.915 g/cm 3 , 
which is the density where the perfect tetrahedral order occurs. The 
behavior of the minimum of D, D m [ n (T), complements the known be- 
havior of D max (T) for the same model [30, 31, 35], where a maximum 
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Figure 3. Dependence of the diffusion constant D on p along isotherms (for p < 
1.0 g/cm 3 ). Open symbols are the new simulations we report, and filled symbols are 
from Ref. [30]. Adapted from [40]. 



occurs due to breaking hydrogen bonds at high pressure; the density of 
the D m i n (T) increases slightly with increasing T, while the density of 
D max (T) decreases with increasing T [31]. This is expected, since the 
range of densities where anomalous behavior occurs expands with de- 
creasing T. We show the loci of D m \ n {T) and D max (T), along with the 
spinodal and locus of density maxima in Fig. 4. 

Below the spinodal, D also increases, since the mobility of the gas is 
larger than that of the liquid. However, the simulations clearly show 
that Anin for the liquid occurs prior to the onset of cavitation, and so 
the location of D ni { n we estimate is not affected by phase separation. Re- 
cently, Ref. [35] estimated the location of D m \ n for the same model along 
several isotherms, and associated D m [ n with a maximum in orientational 
order. 
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Figure 4- Relation of the loci of maxima and minima of D with Tmd and the 
spinodal. Open symbols are from the present work, and filled symbols are from 
Ref. [30]. Adapted from [40]. 

4. Conclusions 

Water exhibits a very complex structure and its properties and anoma- 
lies are strongly influenced by variations of pressure. For high densities 
(p > Pmax), water behaves as a normal liquid and the decrease of D with 
increasing pressure is governed by steric effects. For p min < p < /Vaxj as 
the pressure is decreased, the presence of defects and interstitial water 
decrease, the tetrahedral structure dominates, with stronger hydrogen 
bonds. This process reaches its maximum at p = p min « p icc . Fur- 
ther stretching destabilizes the hydrogen bond network, leading to an 
increase in mobility. The locus of D m \ n roughly tracks the spinodal, not 
surprising since the same breakdown of tetrahedral order that gives rise 
to D mm also facilitates cavitation. 
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